A CONVERGENT APPROXIMATION OF THE PARETO OPTIMAL 

SET FOR FINITE HORIZON MULTIOBJECTIVE OPTIMAL 

CONTROL PROBLEMS (MOC) USING VIABILITY THEORY * 

A. GUIGUE * 

Abstract. The objective of this paper is to provide a convergent numerical approximation of the 
Parcto optimal set for finite-horizon multiobjcctive optimal control problems for which the objective 
space is not necessarily convex. Our approach is based on Viability Theory. We first introduce the 
^sj set-valued return function V and show that the epigraph of V is equal to the viability kernel of a 

properly chosen closed set for a properly chosen dynamics. We then introduce an approximate set- 
valued return function with finite set- values as the solution of a multiobjective dynamic programming 
equation. The epigraph of this approximate set-valued return function is shown to be equal to the 
finite discrete viability kernel resulting from the convergent numerical approximation of the viability 
kernel proposed in [4j [5]. As a result, the epigraph of the approximate set- valued return function 
converges towards the epigraph of V. The approximate set- valued return function finally provides 
the proposed numerical approximation of the Pareto optimal set for every initial time and state. 
Several numerical examples are provided. 
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1. Introduction. Many engineering applications, such as trajectory planning 
for spacecraft [7] and robotic manipulators [13], continuous casting of steel [J3], etc., 
can lead to an optimal control formulation where p objective functions (p > 1) need 
to be optimized simultaneously. For a general optimization problem (GOP) with a 
I vector-valued objective function, the definition of an optimal solution requires the 

J> comparison between elements in the objective space, which is the set of all possible 

values that can be taken by the vector-valued objective function. This comparison 
is generally provided by a binary relation, expressing the preferences of the deci- 
sion maker. In applications, it is common to consider the binary relation defined 
in terms of a pointed convex cone P C R p containing the origin [20j . However, in 



this paper, for simplicity, we will only consider the case P = R^, which yields the 
well-known Pareto optimality. The resolution of (GOP) therefore consists of finding 
the set of Pareto optimal elements in the objective space, or Pareto optimal set. In 
general, this set cannot be obtained analytically and we have to resort to numerical 
approximations. The main objective of this paper is therefore to propose a conver- 
gent numerical approximation of the Pareto optimal set for a general finite-horizon 
multiobjective optimal control problem (MOC). By general, we mean that we do not 
make any convexity assumption on the objective space Y (or more generally, on the 
set Y + R+). Indeed, in the case where the objective space (or Y + R+) is convex, 
simple methods such the weighting method can be used to generate the entire Pareto 
optimal set (Theorem 3.4.4, [17, p. 72]). 

When the objective space is not convex, very few approaches to find the Pareto 
optimal set have been proposed. An important line of research is to use evolutionary 
algorithms [HIH], such as genetic algorithms. Also, very recently, an approach [To] 
inspired from the e-constraint method in nonlinear multiobjective optimization 16, 
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pp. 85-95] has been developed for multiobjective exit-time optimal control problems 
where P — R/j_. In this approach, the n-dimensional state is augmented by p — 1 
dimensions, which yields a new single objective optimal control problem. The Pareto 
optimal set of the original problem can be retrieved by inspecting the values of the 
return function of this new problem. The return function of the new problem is fi- 
nally approximated by solving numerically the corresponding (n+p — l)-dimensional 
Hamilton- Jacobi-Bellman equation using a semi-Lagrangian "marching" method. 

In this paper, instead of an exit-time optimal control problem, wc consider an 
optimal control problem over a finite horizon [0,T]. The proposed approach funda- 
mentally differs from [15 . Instead of augmenting the state space and solving the 
resulting augmented Hamilton-Jacobi-Bellman equation, wc define the set-valued re- 
turn function V(-, •) : [0,T] x R™ — > 2 RP [TTJ fTS] as the set-valued map associating 
with each time t £ [0, T] and state x £ R ra the set of Pareto optimal elements in 
the objective space Y(t, x), where Y(t, x) is the set of all possible values that can be 
taken by the vector-valued objective function for trajectories starting at x at time 
t. Hence, the Pareto optimal set for any time t and state x can be obtained just 
by evaluating V at (i,x). We then derive a convergent approximation of V using 
Viability Theory [HIS]. This approximation is a set-valued map with finite set-values, 
called the approximate set-valued return function. Hence, an approximation of the 
Pareto optimal set V (t, x) can be obtained just by evaluating the approximate set- 
valued return function at (i,x). The advantage of using Viability Theory is that it 
provides a framework that allows to deal with problems with minimal regularity and 
convexity assumptions. Hence, it is expected that the proposed approach could be 
easily extended to more general classes of problems than the one considered in this 
paper, e.g., problems with state constraints, etc.. 

More precisely, the first step in the proposed approach is to show that the epi- 
graph of V, i.e., the graph of the set- valued map V + R+, is equal to the viability 
kernel of a properly chosen closed set for some properly chosen dynamics. The next 
step is to introduce an approximate set-valued return function as the solution of a 
multiobjective dynamic programming equation. The epigraph of this approximate 
set- valued return function is shown to be equal to the finite discrete viability kernel 
resulting from the convergent numerical approximation of the viability kernel pro- 
posed in [H [5]. From there, we easily obtain that the epigraph of the approximate 
set-valued return function converges in the sense of Painleve-Kuratowski towards the 
epigraph of V. The multiobjective dynamic programming equation obtained is very 
similar to the one obtained in (10j . where no proof of convergence was provided. 

This paper is organized as follows. In <|2j we detail the class of multiobjective 
optimal control problems considered and define V. In ^J3j we briefly discuss the con- 
cept of optimality in multiobjective optimization and present several useful properties 
related to Pareto optimal sets. In ^4j we show that the epigraph of V is equal to the 
viability kernel of a properly chosen closed set for a properly chosen dynamics. Fol- 
lowing [H |5] , we then propose in $\ a finite discrete approximation of this viability 
kernel. In f|6j we show that the finite discrete viability kernel resulting from this 
approximation is equal to the epigraph of an approximate set- valued return function, 
defined as the solution of a multiobjective dynamic programming equation. From this 
multiobjective dynamic programming equation, we derive in $7] a numerical algorithm 
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to compute the approximate set-valued return function and therefore the approximate 
Pareto optimal set at (£, x). Some numerical examples [TU], for which the Pareto op- 
timal set is analytically known, are provided in SjH] and some conclusions are finally 
drawn in <J9] 

2. A multiobjective finite-horizon optimal control problem ([3]. ) 

In this paper, we will take for || ■ || in R p and R™ the supremum norm. Let B be 
the closed unit ball. 

Consider the evolution over a fixed finite time interval I — [0,T] (0 < T < oo) of 
an autonomous dynamical system whose n-dimensional state dynamics are given by a 
continuous function f (•, •) : R™ X U — > R™, where the control space U is a nonempty 
compact subset of R m . The function f(-,u) is assumed to be Lipschitz, i.e., some 
K f > obeys 

VueC7, Vx 1)X2 eR n , ||f(x 1 ,u)-f(x 2 ,u)||<ii f ||x 1 -x 2 ||. (2.1) 

We also assume that the function f is uniformly bounded, i.e., some M? > obeys 

Vx £ R", Vu £ U, ||f(x,u)|| < M { . (2.2) 

A control u(-) : I — > U is a bounded, Lebesgue measurable function. The set of 
controls is denoted by U. The continuity of f and the Lipschitz condition (2.1 ) guar- 
antee that, given any t £ I, initial state x £ R™, and control u(-) £ U, the system of 
differential equations governing the dynamical system, 

*(«) = f(x( S ),u( S )), t<s<T, 
x(t) = x, 

has a unique solution, called a trajectory and denoted s — > x(s;i,x, u(-)). Let F be 
the set-valued map defined from R™ to R™ by 



F(x)= (Jf(x,u). 
ueu 

The cost of a trajectory over [t, T], t £ I, is given by a p-dimensional vector 
function J(-, -, •) : I x R n X U -> R p , 

J(t,x,u(.))=/ L(x( S ;x,u(-)),u( S ))d S , (2.4) 

where the p-dimensional vector function L(-, ■) : R n x U — > R p , called the running 
cost, is assumed to be continuous. For simplicity, no terminal cost is included in (2.4). 
We assume that the function L is uniformly bounded, i.e., some Ml > obeys 

Vx £ R n , Vu £ U, ||L(x, u)|| < M L , (2.5) 

and that the function L(-, u) satisfies a Lipschitz condition, i.e., some K^ > obeys 
VueC/, V Xl ,x 2 eR n , ||L(x 1 ,u)-L(x 2 ,u)||<i<'L||x 1 -x 2 ||. (2.6) 

Let L be the set-valued map defined from R™ to R p by 

L(x)= U L (^ U )- 

u£U 
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The objective space Y(t, x) for (MOC) is defined as the set of all possible costs ( |2.4[ ): 

y(t,x) = (j(t,x,u(-)),u(-)ew 



From (2.5), it follows that the set Y(t, x) is bounded (by M\JT), and also that 
Y(t,x) C { — (T — £)Ml1} + R+. However, the set Y(t,x.) is not necessarily closed. 

The set-valued return function V(-, •) : I X R™ — > 2 RP for (MOC) is defined as 
the set-valued map which associates with each time t £ I and initial state x G R™ the 
set of Pareto optimal elements in the objective space Y(t, x), where the definition of 
a Pareto optimal element is postponed to <|3| 

V(t,x)=£(cl(Y(t,x))). (2.7) 

The closure in ( |2.7[ ) is used to guarantee the existence of Pareto optimal elements 
(Proposition 3.5^2]). Hence, Vi € I, Vx G R'\ V(*,x) 7^ 0. 



Remark 2.1. Whenp = 1, (2.7) tafces the form 

V(t,x) = \ inf / L(x(s;t,x,u(-)),u(s)) dsl. 
I "(•)£"./* J 

Hence, V(i,x) = {u(t, x)}, where v(-, ■) is the value function for single objective opti- 
mal control problems [51 119) . 

Finally, as V(t,x) C cl(F(t,x)), we have 

F(i,x)c{-(T-i)M L l} + R^. (2.8) 

The objective of this paper is to find a convergent approximation to the Pareto 
optimal set V(0, Xo) where Xo € R™ is some given initial state. 

3. Multiobjective Optimization. For an optimization problem with ap-dimen- 
sional vector-valued objective function, the definition of an optimal solution requires 
the comparison of any two elements yi,y2 in the objective space, which is the set of 
all possible values that can be taken by the vector-valued objective function. This 
comparison is generally provided by a binary relation, expressing the preferences of 
the decision maker. In applications, it is common to consider the binary relation 
defined in terms of a pointed convex cone P C TV containing the origin 20J . 

Definition 3.1. Let yi,y2 G R p - Then, y 1 ^ y 2 if and only ify 2 eyi+P. 



The binary relation in Definition |3.1| yields the definition of generalized Pareto 
optimality. 

Definition 3.2. Let S be a nonempty subset of R p . An element yi G S 
is said to be a generalized Pareto optimal element of S if and only if there is no 
y 2 S S (y 2 7^ Yi) such that yi G y 2 + P, or equivalently, if and only if there is no y 2 
such that yi G y 2 + P\{0}. The set of generalized Pareto optimal elements of S is 
called the generalized Pareto optimal set and is denoted by £(S,P). When P — R^, 
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the generalized Pareto optimal elements are only referred to as Pareto optimal ele- 
ments, and the set £(S, P) is simply denoted £(S). 

An important role in this paper is played by the external stability or domination 
property [TBI PP- 59-66]. 

Definition 3.3 (External stability). A nonempty subset S of HP is said to be 
externally stable if and only if 

Sc£(S,P) + P. 

An immediate consequence of the external stability property is that S + P = 
£(S, P) + P. When P is closed, a sufficient condition for a nonempty closed set 
S to be externally stable is given in Proposition |3.4| Note that this condition is also 
sufficient to guarantee the existence of generalized Pareto optimal elements. 

PROPOSITION 3.4 (Theorem 3.2.10, [18j p. 62]). Let S be a nonempty closed 
subset of TV. If P is closed and S is P-bounded [HI p. 52], i.e., S + l~l —P = {0}, 
then S is externally stable. 

Corollary 3.5. Let K be a nonempty compact subset of TV . If P is closed, 
then K is externally stable. 

In this paper, for simplicity, we will only consider the case P = T\F + . 

4. Characterization of the set- valued return function. In this section, we 
show that the epigraph of the set-valued return function V, i.e., the graph of the 
set- valued map V + R+ , is equal to the viability kernel of a properly chosen closed 
set for some properly chosen dynamics. 

Define the set-valued maps FL ff from R" to R™ x TV by 
FL CT (x) = co( |J {(f (x, u), <rL(x, u))}] , 

and where co(S') denotes the closure of the convex hull of the set S. Observe that 
the set-valued map FL CT takes convex compact nonempty values. Moreover, FL CT is 
bounded by AfpL = max{Mf,AfL} and Lipschitz with Lipschitz constant -KpL = 
max{K: f ,if L }. Denote FL+(x) = FL +1 (x) and FLr(x) = FL" 1 ^). Define finally 
the expanded set-valued map 4> from R x R™ x R p to R x R" x TV by 

M+ ^_ / {1}xFL-(x) if t<T, , s 

W ' X ' Z) ~ \ [0,1] xco(FL-(x)U{(0,0)}) if t>T. [ > 

It is easy to see that <f> is a Marchaud map (Definition 2.2, p. 184, [4]) bounded by 
max{l, AfpL}. 



Consider now the differential inclusion 

(i(s), x(s), z(s)) G 0(t(s), x(s), z(s)) a.e. s > 0, (4.2) 
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and the closed set H = {(i,x,z), t <= [0,T], x G R n ,z G {-(T - t)M L l} + R^_ . 

Proposition 4.1. TTie epigraph of the set-valued map V is equal to the viability 
kernel ofH for <f>, i.e., 

Epi(F) = Viab^CW), 

where Epi(V) is defined as Epi(V") = Graph(V + R?_). 
Proof. First, we prove the inclusion 

Viab ("H) C Graph(V + R^). 

Take (i,x,z) G Viab^CH). Hence, (t,x,z) G H, or £ G i" and z G {-(T-£)M L 1}+R^ . 

• Assume that t = T. Then, as V(T,x) = {0}, it follows that z e R p + = 
V(T,x)+R$.. 

• Assume that £ G [0, T). Let (£(•), x(-),z(-)) be a solution to ( |4.2| with initial 
condition (£, x, z) which remains in "H. By dehnition of 0, (x(-),z(-)) is a 
solution to the differential inclusion 

(x(s),z(s)) G FL-(x(s)) a.e. sG [0,T-t], 
x(0) = x, 

z(0) = z, 

and £(s) = s + 1. Let s' = s + t. For s' G [t, T], define x'(s') = x(s' - t) and 
z'(s') = z(s' — £). Then, (x'(-), z'(-)) is a solution to the differential inclusion 

(x'(s),z'(s')) G FL-(x'(s')) a.e. s' G [i,T], 
x'(t) = x, 

z'(i) = z, 

By the Relaxation Theorem (Theorem 2.7.2, [19j p. 96]), for e > 0, there 
exists u(-) G U such that 

||x'(.) -x(.;i,x,u(-))|| < e and ||z'(-) - z(-; i, (x, z), u(-))|| < e. 

In particular, we get 

z'(T)<z(T;£,(x,z),u(-)) + el. 

Moreover, as (s',x'(s'), z'(s')) G "H for all s' G [t,T], we have 

zV)e{-(T-s')M+R+> 

or 

z'(T) G R* . 
Hence, 

z(T;t,(x,z),u(-)) = z- J L(x( S ';t,x,u(-)),u( S '))d S '>z'(T)-el>-el, 
or 

z + el > / L(x(V; t, x, u(-)), u(s')) ds', 

which implies that z G c\(Y(t, x)) + R*J_ = V(t, x) + R^_ by external stability. 
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Second, we prove the inclusion 

Graph(U + R$.) C Viab^H). 



Take (i,x,z) e Graph(U + R^). Hence, te [0,T], from (2.8 1, ze {-(T-t)M L l}- 



R^ + . 



Assume that t = T. Then, z e R*. Therefore, (t(-),x(-),z(-)) = ( T , x , z ) is 



a solution to (4.2) viable in "H. 

Assume that ie 0, T). We have z = z' + d, where z' e V(t, x) and d e R+. 

By definition of V(t, x), there exists a sequence u„(-) e W such that 



lim / L(x„(s;i,x, u n (-)),u n (s)) ds = z'. 



(4.3) 



Using the Compactness of Trajectories theorem (Theorem 2.5.3, p~9|. p. 89]) 
and by passing to a subsequence if necessary, we can therefore assume that 
there exists (x(-),z(-)) solution on [t,T] to the differential inclusion 

(x(s'),z(s')) G FL+(x(s')) a.e. s' E [t,T] 

with initial condition (x, 0) such that 

lim ||x„(-;i,x,u n (-))-x(-)|| = 0, 

n— f+oo 

and 



lim ||z„(-;i,(x,0),u„(-))-z(-)|| = 0. 



(4.4) 



From (4.3) and (4.4), we deduce that z(T) = z'. 



Define now (x'(-), z'(-)) as follows: 
(i(s),x'(s),z'(s)) = 



(s + i,x(s + t),z-z(s + i)) s€ [0,T-t], 
(T,x(T),d) s>T-t. 

As z'(T-i) = z-z(T) =z z' = d, it follows that (*(•), x'(-), z'(-)) is a solution of 
( |4.2[ ). It remains to check that (i(s),x'(s),z'(s)) e "H for all s > 0. F or s > T - i, as 
d e R5_, this is obvious. For se [0, T — £], using the definition of z, (4.3), and (4.4), 
we have 

z'(s) = lim / L(x„(s;t,x, u„(-)), u„(s)) ds + d - z n (s + t;i, (x, 0),u n (-)). 



As 



we get 



z n (s + P,t, (x,0),u n (-)) 



.■+/ 



L(x n (s;£,x, u n (-)),u„(s)) ds, 



z'(s)=d+ lim / L(x„(s;£,x, u„(-)),u„(s)) ds. 
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As 



( L(x n (s;t,x,u n (0),u n (s)) ds > -(T- (« + t))M L l 






and d € R>|_, we finally obtain z'(s) > -(T-i(s))M L l. Hence, (£(s),x'(s),z'(s)) e % 

and (£(•), x'(-),z'(-)) is viable in "H. 

D 

Remark 4.1. Proposition |4.1| remains valid if we take for H the closed set 
{(t,x,z), i € [0,T], x e R",z e {-(T - *)M t l} + R^, where M^ > M L . This 
remark will be used in <|6l 

5. Approximation of Viab^'H). In this section, we approximate Viab^H) by 
finite discrete viability kernels. A preliminary step is to approximate Viab^H) by 
discrete viability kernels. Our developments closely follow [4]. 

5.1. Approximation of Viab0("H) by discrete viability kernels. In [4 (The- 
orem 2.14, p. 190), it is shown that Viab^^) can be approximated by discrete viability 
kernels in the sense of Painleve-Kuratowski by considering an approximation <fr e of <j> 
satisfying the following three properties: 

(Ho) 4> e is an upper semicontinuous set- valued map from RxR" x R p to RxR" x R p 

which takes convex compact nonempty values. 
(Hi) 

Graph(0 e ) c Graph(0) + o(e)B where lim g(e) = + . 

e-fO+ 

(H 2 ) 

V(t e ,x £ ,z e ) e RxR™ xR p , |J 0(i,x,z) C <j) e (t e ,x e ,z e ), 

||(t,x,z)-(t e ,x e ,z e )||<Me 

where M = max{l, Mfl} denotes a bound for cj) and e > is the time step discretiza- 
tion. 

Define the set- valued map (f> e , e > 0, from R x R n x R p to R x R™ x R p by 

f{l}x(FL-(x)+eXMB) if t <T - Me, 

cPe{l,,) ~\[0,l]xco{(FL-(x) + eKMB)U{(0,0)}) if t>T-Me. [ > 



Theorem 5.1. The set-valued map <j) e satisfies (Ho), (Hi), and (H2). 
Proof. 
(H ) This follows from the properties of the set-valued map FL~ and the fact that 
FLr(x) + eKMB C co((FL-(x) + eKMB) U {(0, 0)}). 

(Hi) This relation holds with g(e) = emax{l, K}M. 

Let (f e ,x e ,z e ) e R x R" x R p and (s e ,f e ,l e ) € e (i e ,x e ,z e ). Consider the 
following two cases: 
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• t e < T - Me. Then, s e = 1 and (f e ,I e ) e FL _ (x e ) + eKMB. Hence, 
g(e) = eKM. 

• t e >T- Me. We have s £ € [0, 1] and (f £ , l e ) e co((FL~ (x £ ) + eKMB) U 
{(0, 0)}). There exists t > T such that \t-t e \ < Me. We have 

(f 6 ,l e ) G co((FL-(x e ) + eXA/B)U{(0,0)}), 
C co(FL-(x e )U{(0,0)}) + eJ fi:MB. 

Hence, g(e) = emaxjl, K}M. 

(H 2 ) Let (t e ,x e ,z e ) e R x R" x RP. Take (*,x,z) G R x R" x RP such that 
||(t,x,z) — (t e ,x e ,z 6 )|| < Me. In particular, ||x — x £ |j < Me. Consider the 
two following cases: 

• t e < T - Me. This implies t < T. The Lipschitz property of FL~ yields: 
<f>(t,x,z) = {1} x FL-(x) c {1} x (FL-(x e ) + eKMB) = </> £ (£ £ ,x £ ,z £ ). 

• £ £ > T — Me. Then, either t < T. In such a case, as above, 

cf>(t,x,z) = {1} xFL-(x) c {1} x (FL-{x e ) + eKMB) C </> £ (£ £ ,x £ ,z £ ). 

Either t > T. In such a case, using the Lipschitz property of FL~, we 
have 

FL-(x)U{(0,0)} C (FL-(x £ )+e^TAfB)U{(0,0)}. 

Hence, 

0(t,x,z) = [0, 1] x co(FL-(x) U {(0,0)}) c £ (t £ ,x £ ,z £ ). 

D 

5.2. Approximation of Viab^H) by discrete finite viability kernels. In 

[1] (Theorem 2.19, p. 195), it is shown that Viab^(H) can be approximated by finite 
discrete viability kernels in the sense of Painleve-Kuratowski by considering an ap- 
proximation L £ _h of G £ (-) satisfying the following two properties: 



(H 3 ) 



Graph(r £ h ) C Graph(G £ ) + ^(e, h)B where lim ti^l = +. 

e— J-0+,^— 5-0+ £ 



(H 4 ) V(i, l ,x h ,z h )eR,xR5JxRP, 

|J (G £ (i £ ,x £ ,z £ ) + / l B)nR /l xR; i xR p i cr £ ,, l (4,x h ,z h ), 

||(t e ,x e ,z e )-(t h ,x h ,z h )||<' 1 

where h > is the state step discretization, R/i is an integer lattice of R generated 
by segments of length h, G £ is the set-valued map from R x R™ x R p to R x R" x R p 
defined by 

G £ (t £ ,x £ ,z £ ) = {(t £ ,x £ ,z £ )} + e0 £ (t £ ,x £ ,z £ ), 

and r £ ./j is the set-valued map from R^ x RJJ x R^ to R/j x RJJ x R^ defined as follows. 
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• If t h <T-Me - h, r e , h (t h ,x h ,z h ) = 

[t h +e-2/»,t fc +e+2/i]nR fc x({(x h ,z h )}+eFL-(x h )+a £ , fc B)nR2xRj. (5.2) 

• If t h >T- Me -h,T e>h (t h ,ic h) z h ) = 

[i h ,i, l +e+2/i]nR /l xco(({(x h ,z h )}+eFL-(x h )+a e , /l B)U({(x h ,z h )}+2/iB))nR^xRP, 

(5.3) 
where a e ,h = 2h + ehK + e 2 KM. 

We assume that e > 2ft,. 

Theorem 5.2. The set-valued map V Cj f l satisfies (H 3 ) and (H 4 ). 
Proof. 
(H 3 ) This relation holds with ip(e, h) = Ah + ehK, which verifies 

(A + eK)h , 
lim ^^ ^=0+. 

e-*-0+,|->0+ e 



Let (t h ,x h ,z h ) e Rh x R£ x R£ and (s/,,fh,lh) € r e)h (i/,,Xh,z h ). Con- 
sider the following two cases: 

• t h <T - Me-h. Then, s h £ [t h + e - 2h, t h + e + 2h] nR h dt h + e + 
[~2h,2h], and 

(fh,lh) € ({(x h ,z h )} + eFL-(x h )+a 6j/l B)nR£xR£, 
C {(x h ,z h )} + eFL"(x h ) +a £jft B, 
= {(x h , z h )} + e(FL- (x h ) + eKMB) + (2h + ehK)B. 

Hence, (s h J h ,l h ) € (i/ l ,x h ,z h )+e0 e (t/ l ,x h , z h )+max{2/i, 2/i+e/ili'}B = 
G e (ifc,x h ,zh) + (2h + ehK)B. 

• <h > T - Me - h. We have s/, e [t ft , i h + e + 2h] n R^ and (f h , l h ) e 
co(({(x h ,z h )} + eFL-(x h ) + a 6 ,/,B) U ({(x h ,z h )} + /jB)) n R£ x R*. 
There exists t f >T -Me such that < t e - t h < h. We have 

s h e [£ h ,^ + e + 2/i]nR h , 
C [tfc,t h + e + 2ft], 
C £ e + e[0, l] + [-2/i,2/i]. 

Moreover, 

(fh, lh) e co(({(x h , z h )} + eFL- (x h ) + a e , ft B) U ({(x h , z h )} + 2MB)) n R£ x R p h , 

C co(({(x h , z h )} + eFL- (x h ) + a e , h B) U ({(x h , z h )} + 2MB)), 

= {(x h , z h )} + co((eFL-(x h ) + a e , h B) U ({(0, 0)} + 2KB)), 

= {(x h , z h )} + co((e(FL-(x h ) + eKMB) + (2h + ehK)B) U ({(0, 0)} + 2hB)), 

C {(x h , z h )} + co((e(FL- (x h ) + eKMB)) U {(0, 0)}) + {Ah + ehK)B. 

Hence, (s h ,f h ,l h ) G (t f ,x h ,z h )+ecj) e (t < :,x h ,z h )+iaia,x{Ah+ehK,2h}B — 
G e (t e ,x h ,z h ) + (Ah + ehK)B. 
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(H 4 ) Let (t/,,Xh,Zh) € R/> x R£ x R p h . Take (t e ,x e ,z e ) € R x R™ x R p such that 
||(t 6 ,x e ,z e ) — (th,Xh,Zh)|| < ft. In particular, \t e — th\ < h and ||x — Xh|| < ft. 
Consider the two following cases: 

• th < T — Me — ft. This implies t e < T — Me. The Lipschitz property of 
FL~ yields: 

G £ (i £ ,X £ ,Z £ ) + ftB = (i e ,Xe,Z e )+e0e(*6,X e ,Z e )+/lB, 

= (t e + e+[-h,h\) x ({(x £ ,z £ )} + eFL-(x £ ) + e 2 iCA/B + /iB), 

C (t h + e+[-2h,2h]) x ({(x ft ,z ft )}+ftB + eFL-(x/0+eftifB-|-e 2 lO<f + ftB), 

= (t h + e+[-2h,2h])x ({{x h ,z h )} + eFL-(x h ) + a eJl B). 

Hence, 

(G e (t e ,x e ,z e ) + / l B)nR ft x R£ x Rjcr e>ft (4,x h ,z h ). 

• th >T — Me — h. Then, cither t e < T — Me. In such a case, as above, 
G e (t e , x £ , z e )+hB c (t h +e+[-2h, 2ft]) x ({(x h , z A )}+eFL- (x h )+a e , h B). 
As e > 2ft., 

t fc + e + [-2ft, 2ft] c [i ft ,ifc + e + 2ft]. 

Moreover, 

{(x ft ,z /l )}+eFL-(x /l )+ae,/ l BCco(({(x h ,z h )}+eFL-(x h )+a e ,7 l B)U({(x h ,z h )}+2ftB)). 

Hence, 

(G £ (i £ ,x £ ,z £ ) + ftB) nR h x RJJ xR[c r e>/l (^,x h ,z h ). 

Either, i e > T — Me. In which case, using the Lipschitz property of FL~, 
we have 

G e (t e , x £ , z £ ) + frB = (i £ , x e , z e ) + e([0, 1] x co((FL-(x e ) + eJTMB) U {(0, 0)})) + ftB, 

= {[U-h,t t + e + h]) x (co({(x £ ,z e ) + eFL-(x £ ) + e 2 XAfB}U{(x £ ,z e )}) + ftB). 

Now, 

co({(x £ , z e ) + eFL- (x £ ) + e 2 KMB} U {(x £ , z £ )}) + ftB, 

C co({(x h , z h ) + hB + eFL-(x h ) + eKh + e 2 KMB} U ({(x h , z h )} + ftB)) + ftB, 

C co({(x/,, Bfc ) + eFL-(x h ) + a e ,/,B} U ({(x/,, z h )} + 2ftB)). 

Hence, 

(G £ (i £ ,x £ ,z £ ) + ftB)nR ft x RJJ x R p k cr £) 4t Al x h ,z h ). 

D 

6. Convergent approximation of the set-valued return function V. In 

this section, we first introduce a sequence of approximate set-valued return functions 
with finite set-values recursively defined by a multiobjective dynamic programming 
equation [TTJ [T2] . We then show that the epigraphs of these approximate set- valued 
return functions are equal to the sets involved in the calculation of the finite discrete 
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viability kernels of the discrete set Hh = (H + ftB) n R/, x H 1 ^ x R^ for the finite 
discrete dynamics T e f l (Proposition 2.18, p. 195, [4|). This allows us to conclude that 
the sequence of approximate set- valued return functions is finite and that the epigraph 
of the final approximate set-valued return function of this sequence converges in the 
sense of Painleve-Kuratowski towards the epigraph of the set-valued return function 
V. 

Recall the definition oi H = {(*, x, z), t G [0,T], x G R™, z G {-(T - i)M L l} + 
R5_. Here, we take Ml > Ml (Remark 4.1) in the definition of W. Hence, W = 



{(t,x,z), te [0,T], xeR",ze{-(T-t)M L l} + RP. Let 4 = (J + [-M])nRi. 
We define the finite- valued set-valued map V® h from 1^ x R^ to R^ such that 

Graph(^+R^ + )=^, 

where Hh = (H + ftB) n R/j x R^ x R^ . We recursively define now the finite set- valued 
maps V k h , k > 1, from Iy L x R£ to R^ as follows: 

• If t h < T - Me - h, V^ith, x h ) = 

£ ('|(el+a e ^B)nRP+K fc J(</ l +e+[-2^,2^])nR ft ,(x h + e f+a e j l B)nR^), (f , 1) G FL+(x h ) 

• Otherwise, 

V^' 1 (*fc,Xh) = K*fc(*/i,Xh). (6-2) 



(6.1) 



Remark 6.1. The closure in (6.1) is not required anymore as the sets involved 
are finite. 



We aim in the following two propositions to prove that 

fc + 1 +R^ + )cGraph(V^+K / , 



and 

Graph(V/ fc , l + R£ !+ ) = A\ 

where the sets A k are recursively defined (Proposition 2.18, p. 195, |4]) from A = Hh 
and the relation 

A k+l = {(t h , x h , z h ) G A k s.t. T e , h (t h , x h , z h ) niV «}■ 



To simplify the proof of Proposition |6.1| we will assume that T is a multiple of h and 
that e — 2h > 2h, which guarantees that Vt/j e 4, th + e — 2h > ft. We will also 
assume that for all (i/^Xh) € //, x RJJ such that t/i > h, 

V e ° h (t h , Xh ) = {(-(T + ft - t h )M t - h)l}. 

This guarantees that, for all (4,Xh) € I/i X RJJ, for all (f, 1) £ FL + (xh), for all 
4 G (4 + e + [-2ft, 2ft]) n R h , for all x h G (x h + ef + a e>h ) H R£, 

V? h fo, x h ) - {(-(T + ft - 4)M L - ft)l}. 
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We will finally assume that e, h, and Ml have been chosen such that 

eMi + a e ,h<(e-2/i)Mi,. (6.3) 

Proposition 6.1. 

Vfc, Graphic 1 + R£ , ) C Graph(VA + R? +)• 



Proof. We start by proving that this relation holds for k — 0. Take (i/,,Xh, Zh) S 
Graph(V e 1 /l + R^ , ). Hence, Zh £ V^ h {th, Xh) + R^ . . We have two cases to consider: 



1. If t h >T- Me - h. From (6.2|, we directly get z h e V;\(t ft ,x h ) + R£ 



2. Otherwise, t h <T- Me - /T~By (6.1 1, 



V e ] h (t h) x h ) =£^j( e l+a e ^B)nRP+y e h ((i,,+ e +[-2^,2^)nR /l ,(x h +ef+a^ /l B}nR^, (f,l) € FL+(x h )jY 

Let 4 € (4+e+[-2/ l) 2/i])nR h , x h e (jc h +ef+a e ,fcB)nRJ*, z h € V e ° h (f h ,x h ), 
and Zh € (el + a e ,fcB) n R^. Hence, Zh > — (eMx, + a e ^)l. Moreover, from 
above, we have 

zh = {-(T+h-t h )Mi 4 -h)l > (-(T+h-t h -(e-2h)jMj.-h)l = (-(T+h-t h )Mi.-h)l+(e-2h)Mi.l. 

Hence, by ( |6.3[ ), 

zh > (-(r+/»-t h )3J L -ft)l+(e-2ft)3J L l-(cML+a ei /,)l+(cML-l-a e ,fc)l > -(T+^-t, l )I7 L l+(eM L +a e , /t )l. 

Therefore, 

2h + z h € V e ° fc (t A , x h ) + R£ >+ . 

Hence, 

from which we get that 

z h g ^ h (t h ,x h )+R^ >+ c V^ fc (t fc ,x h ) + R^ + +R^ + - V^(t fcj x h ) + R* i+ . 

Assume now that the relation holds up to /c. We aim to prove that it holds for k + 1, 
i.e., 

Graph(T/ fe + 2 + R* >+ ) C Graph^ 1 + R£ + ). 

Take (i h ,x h ,z h ) e Graph(V; fc + 2 + R^ + ). Hence, z h € V^+ 2 (^,x h ) +R* + . We have 
two cases to consider: 



1. If t h >T- Me - h. From (6.2 1, we directly get z h e V*^" (*/,,x h ) + R£ )+ 



2. Otherwise, t h < T - Me - h. By ( |6.1[ ), for some (f , 1) e FL+(x h ), there 
exist f h € (t ft + e + [-2h,2h]) n R h , x h e (x h + ef + a £jh B) n R£, and 
Zh € (el + a ej ftB) n R^ such that 



zhez h + y e \ +1 (f t ,x h ) + R^ 
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From the induction assumption, we get: 

z h ez h + V e k h (t h ,x h )+R p h _ 



z h G |(d+a e , h B)nR^+V e * h ((*h+e+[-2^2ft])nRh 1 (x h +ef+a 6 , h B)nRX) ) (f,l) G FL+(x h )|+R^ + . 



Applying external stability to 



S=i{el+a^ h B)nKl+V e %(t h +e+[-2h,2h])nR h ,( Xh +e{+a eM B)nIll), (f , 1) G FL+(x h )|, 



and using (6.1 1 yields z h € V e J x (t fc , x h ) + R£ 



Proposition 6.2. 



Vfc, Giaph(Vh+R p h ,)=A k 



Proof. This relation holds for k = by definition. Assume now that the relation 
holds up to k. We aim to prove that it holds for k + 1, i.e., 

Graph(V; fc + 1 +R^ + ) = A fe + 1 . 
First, we prove the inclusion 

Graph(^+ 1 +R^ + )cA fc+1 . (6.4) 

Take (t/,,Xh) G IjxRJJ and Zh G V r £ fe / j hl (ih,Xh) + R^ ,. We have two cases to consider: 



1. If t h >T- Me - h, then by (6.2 ), we have (t h , x h , z h ) G Graph(V e fc , l + R£ >+ ). 
Therefore, from the induction assumption, (£/i,Xh,Zh) G A fe . Moreover, 



by (5.3), we also have (t/,,x h ,Zh) G r e]h (t fc ,x h ,z h ). Hence, (^,x h ,z h ) G 



r eih (i^Xh,z h )nA fc ,orr ei , l (i h ,Xh,Zh)nyl' £ ^ 0, which shows that (tfc,x h ,z h ) G 

A fe+1 . 



2. Otherwise, t h < T - Me - h. By fl6T[ ), for some (f , 1) G FL+(x h ), there 
exist t h G (t h + e + [-2h,2h]) n R ft , x h G (x h + ef + a e ,/>B) n R£, and 
Zh G (el + a e ,/iB) H R^ such that 

z h G 2h + V£ fc (t fc , x h ) + R£ >+ . 

Hence, from the induction assumption, (t/ l7 Xh,Zh — Zh) G A k . To get that 
(^,Xh,Zh) G A k+1 , it remains to prove that (th,Xh,Zh) G A k and (£/i,Xh,Zh— 
Zh) G r e> h(t/i,Xh,Zh) where IVfe(fa,x h ,z h ) is given by ([572]). (*h,Xh,Zh) G 



A comes from Proposition 6.1 and the induction assumption, i.e. 



(t h ,x h) z h ) G Gra P h(V; fe + 1 + R£ + ) C Graph(^ + R^ + ) = A k . 
Moreover, we have: 
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(a) t h £ (t h + e+[-2h,2h])nR h , 

(b) S h e(x h + 6f + a fA B)nR^ 

(c) and 

Zh e (el + Oc,fcB) n R£ 

^ -z h e (-el + a e , h B) n R£ 

=*> Zh-Zh G (z h -el + a £ ,fcB)nR£. 



Hence, r ei?l (t^,x h ,z h ) flit ^6 and < 6.4 1 is proved. 

Conversely, we prove the inclusion 

A k+1 C GraphC^ 1 + R* + ). (6.5) 

Take (tft,Xh,Zh) G A +1 . We have two cases to consider: 

1. If t h > T - Me - h. By definition of A k+1 , (t h ,x h ,z h ) € A k . Hence, from 
the induction assumption, we get Zh e V^ /,(£/,, Xh) + R^ , , and from (6.2), 

Zhe^ +1 (4,x h ) + R^. 

2. Otherwise, th < T — Me — h. Then, there exist (i^,Xh,Zh) € A k such that 
(4,x h ,z h )e r £ ^(ift,x h ,z h ), or: 

(a) 4 e (t fc + e + [-2/i, 2ft]) n R^, 

(b) x h e (x h + ef + a e ,fcB) n RJJ, 

(c) z h e (z h - el + a £ , fc B) n K p h =>■ z h e z h + (el + a e , h B) n R£, 

for some (f , 1) G FL~(xh). From the induction assumption, we have Zh G 
V e k h (t h ,x h ) + n p hi+ . Hence, 

z h g (el + a e ,hB) n R£ + V k h (t h , x h ) + R^ + . 
Applying external stability to 

5= |(el+a £A B)nR5:+y e fc /l (^+e+[-2ft,2/ l ])nR, l ,(x h +ef+« £ , ft B)nR," l ), (f, 1) G FL+(x h )|, 



and using JO) yields z h e V^n^,x h ) + R£ >+ . 

D 

Corollary 6.3. The sequence of approximate set-valued return functions V k h 
is finite. 

Proof. This follows from Proposition |6.2| and the fact that the sequence Ak is 
finite (Proposition 2.18, p. 195, @]). 
D 

We denote k(e, h) the last element of this sequence. 



COROLLARY 6.4. The epigraph of the approximate set-valued return function 
V e ff' converges in the sense of Painleve-Kuratowski towards the epigraph of the 
set-valued return function V , i.e., 

r k{e,h) 



Graphs + Itf ) = lim Graph(KT + R? , )• 

e->0+. --J-0+ 
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Proof. From [3] (Theorem 2.19, p. 195), we have 



Viab^ (H ) = lim Viat>r c h (H h ). 

e->0+, -^0+ 



Moreover, from Proposition [41] we have 

Graph(V + R^ ) = Viab^H). 



Finally, from Proposition 6.2 and [4 (Proposition 2.18, p. 195), we have 

Graph(^ £ ' /l) + R£ >+ ) = A k ^ = Vfa&r. ih (« h ). 
The desired result follows. D 



7. A General Numerical Algorithm. In this section, we present a general 
algorithm to determine the approximate set- valued return function V ^ e ' . As shown 

in Proposition 7. ll to find V e h e ' , it is not needed to compute V k h , k = 0, . . . , k(e, h) 
over their entire domain Ih x R?. 



Proposition 7.1. Vfc > 0, \/t h e I h , t h >T- Me-h- k(e - 2h), Vx h e R£, 
V e k +\t h ,x il ) = V e k h (t h ,x il ). (7.1) 



Proof. For k = 0, (|7.1[) directly follows from (6.2 1. Assume now that (7.1| 



holds for k > 0. We prove that (7.1) also holds for k + 1. Let th <= i/u ^ > 



T-Me-h-(k+l)(e- 2h) and */, e (t h + e + [-2h, 2h)) n R?,. Then, 

t h >t h + {e-2h) > T - Me- h- (k + l)(e-2h) + (e-2h) = T-Me-h- k(e-2h). 

From the induction assumption, we get Vxh € RJJ, V(f , 1) € FL+(xh), Vxh € (xh + 
ef + «aB)nR^, 

V k + 1 (t h ,x h ) = V e k h (t h ,5t h ). 
Hence, using (6.1 ), we obtain V e '^ 2 (t/i,Xh) = V e ij" (tfc,Xh), which completes the proof. 



We now present a very general numerical algorithm to approximate the set V(0, xq), 



where for simplicity, we take the initial state Xq in RJJ. From Corollary 6.4 the sug- 



gested approximation to 1^(0, Xo) is given by the finite set V e ^ (— h, Xq)- The pro- 
posed numerical algorithm is composed of two stages. In the first stage (Algorithm 
1), the computational domain is determined using the bound on the dynamics. Once 
the computational domain has been determined, in the second stage (Algorithm 
2), V f h' (—h, xq) is calculated using the multiple dynamic programming equation 



f . h 



( 6.1 )-(|6.2|) together with Proposition 7.1 



Choose for example €i = 1/2* and hi = 1/2 2 *. Let J be the number of discretization 
steps in hi for the interval [—hi,T + hi], i.e., J = (T + 2hi)/hi + 1 (for simplicity, 
we assume that T is a multiple of hi) and let tj — —hi+j*hi. First, we need to 
determine the computational domains Clj, j = 0, . . . , J — 1. 
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Algorithm 1 

Initialization \/tj, —hi < tj < ej — 3hi, £lj — {xo}. Otherwise, Oj = 0. 

Main loop 

1 Set j = 0. 

2 Repeat 

2.1 Set Xh to the first grid point in Clj. 

2.2 Repeat 

2.3 For all tj', t + a - 2h % < ty < tj + a { + 2h { , 

n r = n r U {(x h + e,f + a ei , hi B) n Rl t , (f, 1) G FL+(x h )}. 

2.4 Until all the grid points in flj have been visited. 

3 Until j = J - 1. 

Algorithm 2 

Initialization Vtj, T - Me, - hi < tj < T + hi, Vx h G % V^.'^^Xh) = {0}. 
Otherwise, V e k ^ M) (t^Xh) = 0. 

Let j* be the largest index such that tj- < T — Me, — hi. 

Main loop 

1 Set j = j*. 

2 Repeat 

2.1 Set Xh to the first grid point in flj and A = 0. 

2.2 Repeat 

2.3 For all tj>, tj + e* - 2h % < tj, < tj + a + 2h { , 

A = Au{(ed+a etM B)nKl+V^ e h ' M (t f , (x h +e i f+a 6i , hi B)nRJJ i )), (f,l) e FL+(x h )}. 

2.4 SetV^^Xh)-^). 

2.5 Until all the grid points in flj have been visited. 

3 Until j = 0. 

To reduce the size of the set A in Algorithm 2, it is possible to only keep the 
Pareto optimal elements at each iteration in Step 2.3. This procedure is justified by 
the two following lemmas. 

Lemma 7.2. Let Si and S 2 be two finite subsets o/R p . Then, £(Si U S 2 , P) = 
£(S 1 U£(S 2 ),P). 

Proof. Take zi e £ (Si U S2,P). Assume for contradiction that zi ^ £(Si U 
£(S 2 ,P),P). Then, by external stability, there exists z 2 G £(Si U£(S 2 ,P),P) C 
S 1 ! U S 2 such that zi e z 2 + -P\{0}. But, this contradicts zi e 5 (Si U S 2 ,P). 
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Conversely, take Zi G £ (Si U £(82), P)- Assume for contradiction that Zi ^ £ (Si U 
S 2 ,P)- Then, by external stability, there exists z 2 G £ (Si U S 2 , P) C 5i U S 2 sucri that 
zi € z 2 + P\{0}. Assume that z 2 ^ £2- Then, necessarily z 2 £ Si U £(S2,P)- But, 
this contradicts z± G £(Si U5(5 2 ,P),P). Assume now that z 2 e 5 2 . Then, by exter- 
nal stability, there exists Z3 € f (5 2 , P) such that z 2 G Z3 + P\{0}. Hence, Zi e Z3 + 
P\{0}withz 3 G 5iU£(S , 2 ,P). But, this again contradicts zi G £ (SiU£ (S' 2 ,P),P). D 

Proposition 7.3. Let Si,..., Si be finite subsets of R p . TVien, 

where Ei is recursively defined by Pi = £(Si,P) and the relation 

E t+ i=£(S l+ iUE t ,P). 

Proof. We proceed by induction. For 1=1, this is by definition. Assume now 
that the relation holds up to I. We aim to prove that it holds for / + 1, i.e., 

,1+1 X 

£[\JSij =E I+ i. 

^ i=l ' 



Apply Lemma 7.2 to U»=i &i an d <^i+i- Then, we get 

£({JS u p) = s(si+iUs(\JSi,P),P 



Using the induction assumption, this yields 



■1+1 
£( U^' P ) =£(Si+iUEi,P)=E i+1 . 

i=l 



D 



Using Proposition 7.3 we can change Step 2.3 in Algorithm 2 to Step 2.3' as 

follows: 

2.1 Set Xh to the first point in fij and A = 0. 

2.2 

2.3' For all £,-/, tj + e, - 2h % < t f < tj +e t + 2h u 

A = £(Au{(eM^. M B)nRl+V^;; t M) (t r , (x h +e i f+a ei , ;ij B)nR^)) ! (f,l) G FL+(x h )}). 

2.4SetV e f h M (t J ,K h )=£(A). 

2.5 Until all the points in ftj have been visited. 

8. Numerical Examples. In this section, the algorithms from ^7] are applied 
to a simple class of optimal control problems for which the set-valued return func- 
tion V can be obtained analytically. The convergence of V '. f" (—hi,x ) towards 
U(0,xo) is investigated. Recall that incrementing i by 1 means dividing the time step 
discretization by 2 and the state step discretization by 4. 
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8.1. Description. Consider the following simple autonomous biobjective (p = 
2) optimal control problem. The one-dimensional (n = 1) dynamics is simply 

i(«) = «(*), se[0,T], 

with U = { — 1,1} and initial condition xq. The cost of a trajectory x(-) over / is 
given by 

Ji(0, xq,u(-)) — / P(x(s))u(s) ds and J2(0, xq, "(•)) = / u ( s ) ds, 



JO JO 

where P(«) is a given polynomial. 

The objective space F(0, Xo) for the problem above can be easily determined. Let 
a = /.i({t E I, u(t) = 1})/T where fi(-) denotes the Lebesgue measure, then 

J 2 (0, x a , O) = / u(s) ds = aT - (1 - a)T = (2a - 1)T, 
Jo 

and 

Ji(0,s o ,«(0) - / P(z(s))«(s) ds = / **(*(«))*(*) ds = [0(«(*))]o. 

Jo Jo 



where Q(-) is an antiderivative of P(-)- With 

x(T) = x(0) + / u(s) ds = x + (2a - 1)T, 
Jo 

and defining (5 = (2a — 1)T, we get 

Ji(0, x , u(-)) = Q(x + <5) - Q(x ) and J 2 (0, x , u(-)) = S. 
The set y(0,xo) can finally be obtained by varying 6 between —T and T. 

It is possible to derive a systematic procedure that gives an interval (c [—T,T]) 
such that for all 8 in this interval, the corresponding element (Ji(0, xo, «(•)); -^(O, Xo, u(-))) 
(Q(xo + S) — Q(xo),S) in the objective space is a Pareto optimal element. We will 
not present this procedure here, and therefore assume that the Pareto optimal set 
V(0, Xo) is known. 

8.2. Results. We consider four different polynomials P(-) and initial conditions 
x with T = 0.5, which yield four problems (MOC1), (MOC2), (MOC3), and (MOC4). 
We have chosen these polynomials such that the Pareto optimal sets V(0, x$) present 
different characteristics. For each problem and for i = 3, % = 4, and i = 5, we compute 
V e .fr.' (—hi,xo) using Algorithm 1 and Algorithm 2 from 971 provide the cardi- 
nality of V t ' fr.' (—hi,xo), i.e., \V e . fc" (— hi,xo)\, calculate the Hausdorff distance 

[1 p. 365] H(V^ hi) (-h i ,x ),V(0,Xo)) between V^ hi \-h u X Q ) and V(0,x o ), and 
finally generate a "normalized" Hausdorff distance 

Wv He iM( h v v(n ss n(v^: h '\^,x ),v(o,x )) 

H(V ei ) hi >(-h l ,x Q ),V(0,x Q ))= 

H ( V e 3 M '(-h 3 ,xo),V(0,xo)) 
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Our results are summarized in Tables [8TJ |8.2[ [873| and |8.4| and also in Figures |8T||8.2| 
|8.3[ and |8.4| Tables |8.1[ |8.2[ |8.3| and 8.4 also contain the size of the corresponding 
problem. The first number corresponds to the total number of grid points, i.e., the 
cardinality of all the sets fij . The second number corresponds to the total number 
of successors for all the grid points, where a successor to a grid point Xh is defined 
as a grid point that can be reached from Xh. Using Graph Theory terminology, the 
size of a problem would correspond to the number of nodes and vertices respectively. 
Figures [8T| |8.2| [83} and |8.4| display the objective space, the Pareto optimal set, and 
the approximate Pareto optimal sets for % = 3, i = 4, and i — 5 or each problem. 



(MOC1) P(x) 



1, xq = 1. The set Y(0,Xo) + R/j_ is convex. Hence, it is 



possible to obtain every element of V(0,Xq) using the weighting method. 
(MOC2) P(x) =-x + l, x = 1.5. The set Y(0,x ) - R+ is convex. Only the two 

Pareto optimal elements for 8 = —T and S — T can be obtained using the 

weighting method. 
(MOC3) P{x) = -2x 3 - 15/ 'Ax 2 + 2/Tox + 1/5, x = 0. The set Y(0, x ) + R 2 + is 

nonconvex and the set 1^(0, xo) is nonconnected. More precisely, 1/(0, xo) is 

the union of two sets. 
(MOC4) P(x) = -3/2* - 1/8, x a = 0. This problem is similar to (MOC3). 



Remark 8.1. To reduce the size of the problems, we have proceeded to some 
simplifications in Algorithm 1 and Algorithm 2. First, we have individually com- 
puted a €i ^ f or the dynamics f(-, •) and each component of the running cost L(-, ■). 
Second, we have reduced the interval [tj + Cj — 2hj,tj + £j + 2h%] 
tj + Cj — 2hi- Finally, for any given I, we have reduced the set (cjl 
a single element, i.e., the closest lattice element to ejl, which somehow corresponds to 
setting ot eit hi = for the running cost. Hence, we have set M to max{l,Mf} = 1. 

Table 8.1 
Results for (MOC1). 



to the single time 

Va ti ki B)nR? to 





i=3 


i=4 


i=5 


Size 


(306,5897) 


(1630,65093) 


(10422,856445) 


\V^ h f^(-h u x )\ 


10 


33 


130 


n(v^ e h f i \-h i ,x ),v(o,x )) 


0.091227 


0.046550 


0.022605 


n(v^; h ' ] {-h u x ),v{Q,x )) 


1 


0.5103 


0.2478 



Table 8.2 
Results for (MOC2). 





i=3 


i=4 


i=5 


Size 


(306,10961) 


(1630,132125) 


(10422,1826357) 


\v^; hi \-h h xo)\ 


34 


130 


514 


nv^ e h i : hi \-h i ,x ),v(o,x )) 


0.051067 


0.033192 


0.016627 


H(V e T h i ; hi \-hi,xo),V(0,xo)) 


1 


0.65 


0.3256 
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Table 8.3 
Results for (MOC3) . 





i=3 


i=4 


i=5 


Size 


(306,6529) 


(1630,66613) 


(10422,834285) 


\V^(-h t ,x )\ 


3 


21 


99 


n(v e ^\-h tl x ),v(o,x )) 


0.765685 


0.054420 


0.035360 


n(V e %: M (-hi,x o ),V(0,x o )) 


1 


0.0711 


0.0462 



Table 8.4 
Results for (MOC4). 





i=3 


i=4 


i=5 


Size 


(306,7553) 


(1630,85213) 


(10422,1134221) 


\V^'^{-hi,x )\ 


9 


33 


129 


n(v e k X' hi) (-h i ,x ),v(o,x )) 


0.033857 


0.028646 


0.014031 


n(K%'; hi) (-h i ,x ),v(o,xo)) 


1 


0.8461 


0.4144 




0.15 



Fig. 8.1. (MOC1): Objective space (plain line), Pareto optimal set V(0,a;o) (bold line) and 
approximate Pareto optimal set V ?*' * (— hi,xo) for i = 3 (o), i = 4 (+), and i = 5 (■). 
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Fig. 8.2. (MOC2): Objective space (plain line), Pareto optimal set V(0,a;o) (bold line) and 
approximate Pareto optimal set V ?*' (— h{,xo) for i = 3 (o), i = 4 (+), and i = 5 f-J. 



•Si.'Ji 




0.04 



Fig. 8.3. (MOC3): Objective space (plain line), Pareto optimal set V(0,a;o) (bold lines) and 
approximate Pareto optimal set V f 1 ' (— h{,xo) for i = 3 (o), i = 4 (+), and i = 5 (^-J. 
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-0.5 





- Y(0, ,T ) 


• 


V(0, x Q ) (part 1) 
■ V(0 Tn) fnart 2^ 


o 


4 = 3 


+ 


4 = 4 


• 


4 = 5 



-0.25 



-0.2 



-0.15 



0.05 



Fig. 8.4. (MOC4): Objective space (plain line), Pareto optimal set V(0,:ro) (bold lines) and 
approximate Pareto optimal set V f" * (— hi,xo) for 4 = 3 (o), 4 = 4 (+), and 4 = 5 (■). 
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8.3. Discussion. Because the dynamics and the final time T are the same for 
the four problems, it is normal that the total number of grid points remains the same. 
On the other hand, the running cost is different for the four problems, hence for a 
given grid point Xh, the sets FL + (xh) differ, which explains why the number of suc- 
cessors varies between the four problems. 



The large value, i.e., 0.765685, in Table 8.3 for (MOC3) comes from the fact that 
the approximate Pareto optimal set is not able for i = 3 to capture the upper part of 
the Pareto optimal set. However, as i increases, the approximate Pareto optimal set 



now captures the upper part of the Pareto optimal set, and as result, the Hausdorff 



distance H(V ?' lz '(— hi,x ), V(0,x Q )) decreases considerably. 



As i increases, as expected, for the four problems, a better approximation of the 
Pareto optimal set is obtained. The proposed approach also works well regardless 
whether the set Y(0, xq) + R/j. is convex or not. 

9. Conclusion. In this paper, we have derived a convergent approximation of 
the Pareto optimal set for finite-horizon multiobjective optimal control problems. 
Several techniques such as domain decomposition 3J or dynamic grid refinement @] 
could be considered to improve the computational complexity of the proposed ap- 
proach. Let us also mention the idea of clustering recently proposed in [10], which 
consists of keeping only a subset (with fixed cardinality) of the approximate Pareto 
optimal set at each grid point during the resolution of the multiobjective dynamic 



programming equation (6.1)-(6.2) 



A direct extension to this work would be to consider the general case of a pointed 
closed convex cone P instead of the nonnegative orthant R/L and constraints on the 
state. Also, the proposed approach could be considered for other classes of optimal 
control problems, such as multiobjective exit-time optimal control problems 15j. 
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